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ABSTRACT 

The most widely known modem method for estimating gravity field values from observed data is least- 
squares collocation. Its advantages are that it can make estimates at arbitrary locations based on irregu ar y 
spaced observations, and that it makes use of statistical information about errors in the input data while pro- 
viding corresponding information about the quality of the output estimates. Disadvantages of collocation in- 
clude the necessity of inverting square matrices of dimension equal to the number of data values and the need 
to assume covariance models for the gravity field and the data errors. Fourier methods are an important alter- 
native to collocation. They have the advantage of greater computational efficiency, but require data estimates 
to be on a regular grid and do not use or provide statistical accuracy information. 

The GEOFAST algorithm is an implementation of collocation that achieves high computational effi- 
ciency by transforming the estimation equations into the frequency domain where an accurate approximation 
may be made to reduce the workload. The forward and inverse Fast Fourier Transforms (FFT) are utilized. 
We have demonstrated the accuracy and computational efficiency of the GEOFAST algorithm using two sets 
of synthetic gravity data: marine gravity for an ocean trench region including wavelengths longer than 
200 km; and local land gravity containing wavelengths as short as 5 km. We discuss these results along with 
issues such as the advantages of first removing reference field models before carrying out the estimation a go- 

rithm. 


1. INTRODUCTION 

Algorithms for estimating gravity field quantities should be theoretically sound, accurate, and computa- 
tionally efficient. In characterizing and describing such algorithms, the issues we consider mclude the specific 
geodetic quantity or quantities to be computed or estimated, the general approach (deterministic or statisti- 
cal), the type or types of measurement data to be used, the geographic distribution of points at which geodetic 
quantities are observed and are to be computed, and the approximations (if any) made in the computation 
algorithm. The class of least-squares collocation methods is statistical in nature; general m terms ol the quan- 
tities to be estimated, the measured data to be used, and the geographical distribution of data and computation 
points; and, in principle, exact if assumed covariance models are reasonable. 

GEOFAST is an approach to implementing collocation with high efficiency by carrying out computa- 
tions in the spatial frequency domain. It assumes the availability of two-dimensional gndded data; it uses 
stationary statistical gravity field models; it computes minimum-variance estimates of gravity field quanti- 
ties; and it achieves high computational efficiency of the order of N log N operations, where is t e num 
of measurement points. This efficiency depends on the use of a critical approximation which, in practice 
introduces little loss of accuracy. Furthermore, the user has control of the tradeoff between accuracy and 
computational efficiency. 

2. THE GEOFAST ALGORITHM 

The GEOFAST algorithm provides an efficient computational solution to the minimum-variance estima- 
tion equations of collocation: 

x = c«(c zz r 1 z (1) 
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, ,. . xj r , nf Himemion N bv N) is the sum of its autocovariance plus 

where z is a data vector of dimension N, C zz (ol dimension in y > 

the order of N 3 + MN operations, limiting the feasible applications of collocauon methods where data 


on 


of thousands of measurements are involved. 

This workload can be reduced to the order of N log N, if the following assumptions are valid: measured 

GEOFAST algorithm. 
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Figure 1 . Outline of the GEOFAST algorithm. 

Transformation of equation (1) to the spatial frequency domain is carried 
of the Fast Fourier Transform (FFT), leading to an equivalent estimation equation in 

ables 2L' and L': 

x' = C' xz (C' zz ) _1 z' 

practice, the firs, step of the GEOFAST algorithms 

lions^how^ver ^e^t eaneL^Ppmxhnations^sMtrniaUo the computational economy of the GEOFAST algo- 
rithm are"incorporated. In effecMhe Cn' ami C XI ' matrices that would be the exact transforms are 
banded approx, ma.ions, and a compensating data windowing is applied tn the 
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computation of x ' . The final step (right-hand column of Figure 1 ) is the inverse FFT from x ', the frequency- 
domain solution, to x, the spatial domain solution, incorporating the inverse of the windowing operation that 

was part of the transformation from ^ to z' . 

User control of tradeoffs between accuracy and computing time in the GEOFAST algorithm involves the 
selection of parameters governing the bandwidth retained in the approximations to the matrices C zz ' and C xz ', 
and the number of iterations carried out in the iterative solution of the approximate frequency-domain equa- 
tions. In actual applications, reasonable parameter choices result in an accurate and highly efficient algorithm. 

3. FIRST TEST DATASET — MARINE GRAVITY 

For the first set of tests, gridded values of gravity anomaly and deflection of the vertical were generated 
for a 36 deg by 36 deg area in the western Pacific that includes a significant ocean trench, using the Rapp 
(1981) worldwide spherical harmonic expansion to degree and order 180. The GEOFAST algorithm is used 
to compute components of the deflection from the gravity anomaly data; comparison with the deflection com- 
ponents computed originally from the spherical harmonic model (regarded as truth data) provides a measure 
of algorithm performance. To quantify the effects of using high-order reference fields, the same spherical 
harmonic model can be used to generate a set of anomalies and deflections, but for degrees 91 through 180 
only. Using GEOFAST with this high-frequency field is equivalent to removing a degree-and-order-90 refer- 
ence field from the data prior to processing. Other parameters explored in the testing include data extent, grid 
spacing, GEOFAST bandwidth, choice of covariance models, and measurement noise level. 



mo ” o 


Figure 2. Input data and results for GEOFAST test. Input gravity anomalies are shown at the 
left; synthetic deflection truth data are shown in the center; errors in the estimated deflection 
appear on the right. 

A typical test is illustrated in Figure 2. The gravity anomaly data (in this case with the degree-90 reference 
field removed) at 0.5 deg grid spacing are shown on the left; these are used as input to GEOFAST. Deflection 
data computed by GEOFAST for the same grid, assuming an attenuated white noise covariance model (with 
correlation distance of 0.5 deg), will be compared to the synthetic “truth data” shown in the middle. The re- 
sulting deflection error field appears on the right. Some test results are summarized in Table 1 , which shows 
the effects of the bandwidth parameter, m B , and the grid spacing (0.25 or 0.50 deg) on accuracy and comput- 
ing time. Table 2 summarizes the effects of measurement noise on estimation accuracy. To a first approxima- 
tion, measurement noise (in this case, white noise) contributes additively to estimation error variance. T e 
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Table 1. Summary of GEOFAST Test Results. 


1 Resolution 

Bandwidth 

Error (arc* 

>ec) 

Grid (deg) 

Points 

MB 

Time* (min) 

Mean 

Sigma 

Maximum 

0.25 

3600 

0 

4 

0.001 

0.15 

0.70 

0.25 

3600 

4 

30 

0.003 

0.09 

0.31 

0.50 

900 

0 

I 

0.003 

0.17 

0.59 


Table 2. Effects of Measurement Error on GEOFAST Performance. 


Grid (deg) 

Measurement Error 
rms (mgal) 

Noise-Induced Error 
rms (arcsec) 

Total Estimation Error 
rms (arcsec) 

0.25 

0 

0.00 

0.15 

1 

0.03 

0.15 


2 

0.07 

0.17 


5 

0.16 

0.22 

0.50 

0 

0.00 

0.17 

1 

0.04 

0.17 


2 

0.08 

0.19 


5 

0.20 

0.26 


effects of removing a reference field from the measured data were explored in another set of tests. A typical 
result is that the error standard deviation is reduced by a factor of more than five by removing the reference 
field. 


4. SECOND TEST DATASET — LAND GRAVITY 

The second set of tests was designed to investigate the performance of GEOFAST on local data sets 
containing high-frequency information. Synthetic data were generated, using amass dipole model originally 
formulated for the gravity gradiometer testing program. The statistical properties of the data are consistent 
with those for an actual test area in north Texas. As in the marine gravity tests, synthetic truth data were 
generated in the form of gravity anomalies and deflections of the vertical, for each point of a 5 km grid cover- 
ing a square 500 km on a side. GEOFAST is used to estimate the east component of the deflection from the 
anomaly data; comparison with the directly generated deflection values provides a measure of the algorithm s 
performance. These tests use a covariance model with a correlation distance of 10 km. Estimation regions 
ranging from 40 to 150 km on a side are selected from the interior of a 150 km square located centrally within 
the original 500 km area, with resulting error standard deviations between 0.4 and 1 .9 arcsec. The estimation 
accuracy would be improved by using a finer grid spacing. 

5. SUMMARY 

The GEOFAST algorithm, described briefly above, is used to estimate deflection of the vertical from 
gravity anomaly to an rms accuracy of better than 0.2 arcsec with modest computer cost. Conventional imple- 
mentations of collocation require considerably more computing effort. GEOFAST provides a user tradeoff 
between computing cost and estimation accuracy. In a typical application reducing the accuracy requirement 
from 0.1 to 0.2 arcsec reduces computing time by a factor of eight. Other test results confirm the importance 
of removing a high-order reference field from the data before applying the algorithm. For example, removing 
a degree-90 reference field from the synthetic marine gravity data reduces estimation errors by a factor of 

five. 

The GEOFAST algorithm can be applied to current real-world problems involving the reduction of grad- 
iometer or altimeter data. Potential application areas include the production of accurate maps for marine crus- 
tal and lithospheric studies. 
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